library('rjags')
library('coda')
library('jagshelper')
library('jagsUI')
library('mcmcplots')  # caterplot
library('dplyr')

criticism_results

Random effects model

# Model

model_code <- "model
{
  for(i in 1:32){
    r[i] ~ dbin(p[i], n[i])                                     #Likelihood     
    #Log-odds for control (t[i]=1) and treatment (t[i]=2) groups        
    logit(p[i]) <- mu[s[i]] + delta[s[i]]*equals(t[i],2)        
    #Deviance contribution
    rhat[i] <- p[i] * n[i]                                          #  expected value of the numerators 
    dev[i] <- 2 * (r[i] * (log(r[i])-log(rhat[i]))  +  (n[i]-r[i]) * (log(n[i]-r[i]) - log(n[i]-rhat[i])))
  }
  
  for (j in 1: 16){
    #Priors for baseline effects           
    mu[j] ~ dnorm(0,1.0E-6) 
    #Hierarchical random effects model for treatment effects
    delta[j] ~ dnorm(d, prec)       
  }
  
  or<-exp(d)                                    #Population odds ratio
  
  d ~ dnorm(0.0,1.0E-6)             #Prior for population treatment effect
  tau ~ dnorm(0.0,1.0E-6)I(0,)      #Prior for between studies sd
  prec<- 1/(tau*tau)
  
  delta.new~dnorm(d,prec)               #Replicate log OR for prediction
  
  # For plotting purposes put population treatment effect & predicted effect from new study
  # in elements 19 & 20 respectively
  delta[19] <- d
  delta[20] <- delta.new
  
  #Total Deviance
  resdev <- sum(dev[])
}
" %>% textConnection


# Data

data <- "s  t   r   n
1   1   2   36
1   2   1   40
2   1   23  135
2   2   9   135
3   1   7   200
3   2   2   200
4   1   1   46
4   2   1   48
5   1   8   148
5   2   10  150
6   1   9   56
6   2   1   59
7   1   3   23
7   2   1   25
8   1   1   21
8   2   0   22
9   1   11  75
9   2   6   76
10  1   7   27
10  2   1   27
11  1   12  80
11  2   2   89
12  1   13  33
12  2   5   23
13  1   8   122
13  2   4   130
14  1   118 1157
14  2   90  1159
15  1   17  108
15  2   4   107
16  1   2103    29039
16  2   2216    29011" %>% read.delim(text=.) %>% as.list


# Starting/initial values

initial_values <- list(
list(d = 0,  tau=1,mu = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
        delta = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,NA,NA,NA,NA),
        delta.new=0
        ),
        
list(d = -1, tau=2, mu = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
delta = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  1, 1,NA,NA,NA,NA),
        delta.new=0
        ),
        
list(d = 1, tau=0.5,mu = c(-1,-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1),
delta = c(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,  -1, -1, -1, -1, -1,NA, NA,NA,NA),
        delta.new=0
        )
        
)

parameters_to_save <- c("d","delta", "dev", "or", "prec", "resdev", "tau")


results <- jags(data = data,
            inits = initial_values,
            parameters.to.save = parameters_to_save,
            model.file = model_code,
            n.chains = length(initial_values),
            n.adapt = 100,
            n.iter = 50000,
            n.burnin = 20000,
            n.thin = 2)
## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 32
##    Unobserved stochastic nodes: 35
##    Total graph size: 668
## 
## Initializing model
## 
## Adaptive phase, 100 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 20000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 30000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.
summary(results)
##                   mean         sd          2.5%          25%          50%
## d          -0.90404672 0.27436264 -1.498971e+00  -1.07040605  -0.88628981
## delta[1]   -0.89421922 0.66465245 -2.287918e+00  -1.29759367  -0.86638632
## delta[2]   -1.02165360 0.36495353 -1.774609e+00  -1.26019716  -1.00735616
## delta[3]   -1.09008801 0.56090320 -2.294945e+00  -1.44025620  -1.05673516
## delta[4]   -0.73808237 0.69647378 -2.141689e+00  -1.16567311  -0.73457252
## delta[5]   -0.16606609 0.42365551 -9.751880e-01  -0.45184314  -0.17458763
## delta[6]   -1.49490200 0.62043872 -2.861497e+00  -1.86512747  -1.44036861
## delta[7]   -1.02280997 0.65378130 -2.418614e+00  -1.42283670  -0.98497706
## delta[8]   -1.07187369 0.77318832 -2.769949e+00  -1.51865204  -1.01682417
## delta[9]   -0.77681697 0.43550751 -1.659283e+00  -1.06202548  -0.76817809
## delta[10]  -1.37371519 0.63799404 -2.791584e+00  -1.75007628  -1.31932819
## delta[11]  -1.47272171 0.55423410 -2.671270e+00  -1.81654983  -1.43242253
## delta[12]  -0.88374776 0.47939189 -1.867594e+00  -1.19062217  -0.86933493
## delta[13]  -0.85593265 0.48421207 -1.854932e+00  -1.16352476  -0.84337894
## delta[14]  -0.33014205 0.14288969 -6.113890e-01  -0.42544151  -0.32966958
## delta[15]  -1.32469315 0.46513230 -2.321688e+00  -1.61673391  -1.29823410
## delta[16]   0.05515721 0.03156944 -7.050354e-03   0.03395599   0.05534603
## delta[19]  -0.90404672 0.27436264 -1.498971e+00  -1.07040605  -0.88628981
## delta[20]  -0.90453094 0.81484732 -2.603124e+00  -1.37482915  -0.87563495
## dev[1]      0.78735352 1.10807900  7.468822e-04   0.08096559   0.36077211
## dev[2]      0.53477621 0.74678142  5.076015e-04   0.05533416   0.24874425
## dev[3]      0.94435397 1.35321453  8.452687e-04   0.09403138   0.42569557
## dev[4]      0.83235452 1.17899333  8.535195e-04   0.08569379   0.37932149
## dev[5]      0.95839621 1.32632842  9.722126e-04   0.09718985   0.44133241
## dev[6]      0.62435730 0.90658331  5.857951e-04   0.06104878   0.27529279
## dev[7]      0.67774845 0.97127379  6.724926e-04   0.06838231   0.30640112
## dev[8]      0.83807158 1.05226980  1.056953e-03   0.10464266   0.44580652
## dev[9]      1.05947697 1.46884244  1.032915e-03   0.10908266   0.49171403
## dev[10]     1.32284465 1.72753609  1.432552e-03   0.15064627   0.66926838
## dev[11]     1.25017494 1.68448128  1.276456e-03   0.13541751   0.60026459
## dev[12]     1.24572747 1.49305217  1.828778e-03   0.18141763   0.71496885
## dev[13]     0.88625192 1.23092127  9.562462e-04   0.09319689   0.40702141
## dev[14]     0.50900376 0.74751540  5.436790e-04   0.05102252   0.22750336
## dev[15]     1.30184239 1.72533552  1.370873e-03   0.14451621   0.63558217
## dev[16]     0.59705941 0.72151153  1.051376e-02   0.13253888   0.34963979
## dev[17]     0.86366925 1.22443981  8.345482e-04   0.08905054   0.39617209
## dev[18]     0.80133365 1.12969110  7.677557e-04   0.08175295   0.36657847
## dev[19]     1.17542425 1.59741712  1.165849e-03   0.12716474   0.55844176
## dev[20]     1.00506404 1.29062978  1.283094e-03   0.12342378   0.52457885
## dev[21]     1.17367299 1.60919120  1.208699e-03   0.12184268   0.54465025
## dev[22]     1.08040497 1.41836296  1.304302e-03   0.12302133   0.53539430
## dev[23]     0.87511817 1.23982980  8.825754e-04   0.08751937   0.39767095
## dev[24]     0.72395802 1.03681573  7.501261e-04   0.07405602   0.32310103
## dev[25]     0.86180987 1.20707065  8.809494e-04   0.08892404   0.39873539
## dev[26]     0.73426610 1.03158775  7.017103e-04   0.07463550   0.33391356
## dev[27]     0.98315903 1.37606510  9.880295e-04   0.09814564   0.44760834
## dev[28]     1.00101638 1.40784023  1.045017e-03   0.10143130   0.46040310
## dev[29]     1.05414925 1.47132647  1.014276e-03   0.10762229   0.48830996
## dev[30]     0.86773712 1.21393355  8.386151e-04   0.08883158   0.40525211
## dev[31]     0.99740841 1.40256148  8.443172e-04   0.09985963   0.45649128
## dev[32]     1.00157199 1.43044302  9.458723e-04   0.10067012   0.45143105
## or          0.41994720 0.11237226  2.233598e-01   0.34286927   0.41218219
## prec        2.69476693 2.08145531  5.770277e-01   1.36824353   2.13544447
## resdev     29.56955673 7.36229928  1.691502e+01  24.35465133  28.97018614
## tau         0.72180634 0.24855531  3.507458e-01   0.54656991   0.68431468
## deviance  147.35066741 7.36229928  1.346961e+02 142.13576201 146.75129682
##                    75%        97.5%     Rhat n.eff overlap0         f
## d          -0.72026194  -0.41369940 1.000367  6562        0 0.9996000
## delta[1]   -0.46647224   0.36237499 1.000168 22433        1 0.9224444
## delta[2]   -0.77127144  -0.33524492 1.000059 27851        0 0.9984222
## delta[3]   -0.70827617  -0.06770754 1.000092 23757        0 0.9817111
## delta[4]   -0.30514920   0.66051211 1.000119 28083        1 0.8702444
## delta[5]    0.10872347   0.69350707 1.000096 21499        1 0.6610667
## delta[6]   -1.06140024  -0.43354681 1.000136 13541        0 0.9978222
## delta[7]   -0.58746164   0.18071606 1.000092 20868        1 0.9530889
## delta[8]   -0.57100878   0.33369052 1.000516  4712        1 0.9351556
## delta[9]   -0.48654131   0.06232654 1.000087 24886        1 0.9651778
## delta[10]  -0.93465273  -0.27016972 1.000175 15661        0 0.9931778
## delta[11]  -1.08517635  -0.50605425 1.000231  8485        0 0.9993556
## delta[12]  -0.56537799   0.03262513 1.000099 18559        1 0.9706889
## delta[13]  -0.53276981   0.06231826 1.000087 22341        1 0.9656667
## delta[14]  -0.23406314  -0.04927802 1.000049 27370        0 0.9897111
## delta[15]  -1.00343298  -0.48242332 1.000189 11814        0 0.9994000
## delta[16]   0.07633159   0.11704247 1.000173 16750        1 0.9588222
## delta[19]  -0.72026194  -0.41369940 1.000367  6562        0 0.9996000
## delta[20]  -0.41118986   0.67863958 1.000122 13479        1 0.8905778
## dev[1]      1.04026749   3.96167704 1.000611 17062        0 1.0000000
## dev[2]      0.71740205   2.61874164 1.000102 45000        0 1.0000000
## dev[3]      1.23792065   4.77623317 1.000071 45000        0 1.0000000
## dev[4]      1.10038789   4.19497627 1.000507 13008        0 1.0000000
## dev[5]      1.28986786   4.74213813 1.000547  7444        0 1.0000000
## dev[6]      0.81469641   3.17465660 1.000175 45000        0 1.0000000
## dev[7]      0.89285030   3.40195592 1.000519 45000        0 1.0000000
## dev[8]      1.18375576   3.77322027 1.000084 45000        0 1.0000000
## dev[9]      1.41699205   5.22151535 1.000591 13225        0 1.0000000
## dev[10]     1.84536534   6.16032135 1.000074 45000        0 1.0000000
## dev[11]     1.69004286   5.99292284 1.000188 12600        0 1.0000000
## dev[12]     1.77348433   5.38355679 1.000277 15380        0 1.0000000
## dev[13]     1.17966811   4.43787036 1.000463 29021        0 1.0000000
## dev[14]     0.65633988   2.62732306 1.000064 45000        0 1.0000000
## dev[15]     1.78535184   6.16953446 1.000023 45000        0 1.0000000
## dev[16]     0.78553457   2.60474052 1.000728 11160        0 1.0000000
## dev[17]     1.12871020   4.33454997 1.000350 12311        0 1.0000000
## dev[18]     1.06478149   3.97434027 1.000044 45000        0 1.0000000
## dev[19]     1.58625645   5.70697982 1.000181 35523        0 1.0000000
## dev[20]     1.40510905   4.58963068 1.000010 45000        0 1.0000000
## dev[21]     1.59432758   5.78060127 1.000100 45000        0 1.0000000
## dev[22]     1.49878815   5.02284087 1.000160 17824        0 1.0000000
## dev[23]     1.16384069   4.39252972 1.000457 16928        0 1.0000000
## dev[24]     0.94399265   3.71613117 1.000070 45000        0 1.0000000
## dev[25]     1.14929951   4.23245015 1.000613 17251        0 1.0000000
## dev[26]     0.98274602   3.69896492 1.000005 45000        0 1.0000000
## dev[27]     1.31463849   4.94617287 1.000023 45000        0 1.0000000
## dev[28]     1.31390095   5.03227690 1.000184 24879        0 1.0000000
## dev[29]     1.39679081   5.25197603 1.000063 45000        0 1.0000000
## dev[30]     1.15315016   4.28566326 1.000871 11616        0 1.0000000
## dev[31]     1.33160686   4.95317272 1.000385 29059        0 1.0000000
## dev[32]     1.32489573   5.12587001 1.000107 45000        0 1.0000000
## or          0.48662477   0.66119968 1.000227  8105        0 1.0000000
## prec        3.34740738   8.12858564 1.000797  8125        0 1.0000000
## resdev     34.13914734  45.68301988 1.000171 11073        0 1.0000000
## tau         0.85490587   1.31644180 1.000415  7104        0 1.0000000
## deviance  151.92025802 163.46413056 1.000171 11073        0 1.0000000
plot(results)